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ABSTRACT 

A numerical scheme is described for including radiation in multi-dimensional general- 
relativistic conservative fluid dynamics codes. In this method, a covariant form of the Ml 
closure scheme is used to close the radiation moments, and the radiative source terms are 
treated semi-implicitly in order to handle both optically thin and optically thick regimes. The 
scheme has been implemented in a conservative general relativistic radiation hydrodynamics 
code K0RAL. The robustness of the code is demonstrated on a number of test problems, includ- 
ing radiative relativistic shock tubes, static radiation pressure supported atmosphere, shadows, 
beams of light in curved spacetime, and radiative Bondi accretion. The advantages of Ml clo- 
sure relative to other approaches such as Eddington closure and flux-limited diffusion are 
discussed, and its limitations are also highlighted. 
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1 INTRODUCTION 

Accretion disks are the power source behind many astrophysical 
systems. They span a wide range of central object mass, size and 
type, e.g., young stellar objects, neutron star binaries, black hole 
binaries, gamma-ray bursts, active galactic nuclei, to name a few, 
and exhibit a variety of different regimes of physics. In many sys- 
tems, radiation is so intense that it strongly couples to the accret- 
ing gas and dramatically alters the flow structure and dynamics. To 
correctly infer the physics of such systems, we need models that 
properly take into account the interaction of gas and radiation. 

Proper treatment of radiation is especially important for black 
holes (BHs) near the Eddington luminosity limit, Lrad- Transient 
black hole binaries (BHBs) approach near-Eddington accretion 
rates near the peak of their outbursts (IMc Clintock & Re millardl 
2006): iRemillard & McClintockl 12006) : iDone. Gierlinski & Kubot3 



2007), while some exceptional BHB s spend extended periods of 
time with L > L Edd ( e.g., SS433, lMargpnet"al]|l979l ; iMargonl 
1 1984 GRS1915+105, iFender & Belloni |2004|) . Ultra-luminous 
X-ray sources have even larger luminosities and may conceiv- 
ably be highly super-Eddington stellar-mass BHs (Wata rai et al.l 
12001 ). or intermediate mass BHs accreting at close to Edding- 
ton dMiller & Co lbert 2004). In either case, radiation must play an 
important role. Finally, luminous active galactic nuclei, especially 
those whose supermassive BHs (SMBHs) are growing rapidly in 
mass, may be pe rfect examples of systems with closely coupled 
gas and radiation dCollin et al.ll2002l) . 
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Since most of the radiation from an accretion disk originates 
in the inner regions, general relativistic (GR) effects play an im- 
portant role in determining the emergent radiation spectrum. Treat- 
ing the interaction of fluid dynamics and radiation quantitatively in 
GR numerical codes is a daunting task that has been achieved only 
recently and only for optically thick flows, as we discuss below. 
However, the emergent spectrum is established at an optical depth 
of order unity, t ~ 1, which requires a proper treatment of the 
optically-thick disk (t » 1), the optically-thin corona (t <k 1) if 
one is present, and the optically thick-to-thin transition in between. 
In this paper, we present a GR numerical radiation hydrodynam- 
ics technique and a code which handles all three regimes of optical 
depth. 

Depending on the mass accretion rate, a given accretion sys- 
tem can switch between different spectral states, with different 
radiation mechanisms dominating and with varying degrees of 
coupling between radiation and gas. At very low accretion rates, 
L/Lsdd <K 10~ 2 , e.g., Sgr A* the SMBH at our Galactic Center 
dNaravan. Yi. & M ahadevai j 1995b . the accretion flow adverts most 
of the viscously released energy into the BH rather than radiat- 
ing the energy (some energy probably also flows out in a wind). 
Such flows are optically thin, radiatively inefficient and geomet- 
rically thick, and are called advection-dominated accretion flow s 
(ADAFs, see lNaravan & Yilll994fl995l ; lAbramowicz et al.|[l995h . 
Analytical and semi-analytical models are reasonably successful in 
accounting for the main features in the spectra of ADAFs (e.g., 
lYuan. Ouataert. & Naravanl 2003). However, it was realized early 
on that numerical simulations are necessary to understand properly 
the physics of ADAFs. 

The low optical depth and the extreme radiative ineffi- 
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ciency of an ADAF are a major advantage for simulations 
since they allow one to solve the fluid dynamics separately 
from the radiation fielcQ A number of sophisticated GR mag- 
netohydrodynamic (GRMHD) codes have been developed (e.g. , 
De Villiers etai]|2003l : iGammie etal]|2003l ; I Aminos et all 120051 ; 
Del Zanna et al J 12007). and ADAF-like models have been simu- 



lated using these (see Narayan et alj|2012l : McKinnev et alj|2012l : 
iTchekhovskov & McKinneM 20121 for recent work). To compute 
observables, the output from such pure GRMHD codes are usually 
post-processed with stand-alone radiation trans fer schemes (e.g. 
ISchnittman et al.ll2006l : IShcherbakov et al.ll2010h . Alternatively, a 
simpl e local cooling prescription is inc luded directly into the code 
(e.g.. iFragile & Meier||2009l : iDibi et alj|2012h ; this is not difficult 
since the gas is optically thin, so no radiative transfer is involved. 

For higher accretion rates, 10~ 2 < L/L EM < 0.3, efficient 
cooling sets in, and the inner accretion disk colla pses into an 
optically thick geomet r ically thin accretion disk ( Narayan & Yil 
19951: lEsin et all I Il997l fl998l; iMever-Hofmeister & Meverl 120031 ; 
McClintock & RemillardfeOOfjl : iDone. Gierlifiski & Kubotall 20071) . 
This state of accretion is the best understood of all states, and its 
thermal black-body-like spectrum is well-described by the stan- 
dard thin disk model with a viscosity dShakura & Sunvaevll 19731 : 
lNovikov&Thornelll973l : lFranketal.ll2002l) . However, despite its 
many successes, the a disk model cannot describe important micro- 
physical aspects of the flow, e.g., thermal and viscous stability, ver- 
tical structure of the disk, and the role of magnetic fields. Moreover, 
optically thick, geometrically thin disks are much more difficult to 
simulate numerically because of the strong coupling between gas 
and radiation. 

Only a few time-dependent GRMHD numerical simulations 
of thin disks have been computed so far, and all are based on an ad 
hoc (although physically motivated) local cooling function which 
artificially removes excess heat to keep the disk geometri cally thin 



dShafee et alj|2008l: 



Reynolds & Fabian 



Penna ct 



at to keep the disk geometrically thin 
ai] |2010l : lNo"ble et al.ll2oTTI : see also 



20081) . Although this approach has produced 



valuable results on the dynamics of the gas, it misses real physics 
involving propagation of photons along curved geodesies, radiation 
pressure support, radiative winds, etc. 

Using a flux-limited diffusion approximation, small patches 
of radiatively efficient thin accretion disks h ave been simulated 
using the local shearing box approximation ([T urner et al. 2003; 
iKrolik et alj|2007l ; iBlaes et alj|2007l l201ll : iHirose et alj|2009al lbh. 
and even a few full disk simulations have be en done with a non- 
relativistic code using flux -limited diffusion dOhsuga et alj|2009l ; 
lOhsuga & Mineshigell201 lh . However, to model a radiatively effi- 
cient disk self-consistently it is necessary to handle the radiation 
field in both the optically thick (disk interior) and optically thin 
(corona) limits. This is difficult with the flux-limited diffusion ap- 
proximation which artificially enforces the flux to follow the local 
gradient of the radiative energy density. Fueled by advances in ra- 
diation algorithms, more advanced radiation moment closures are 
now becoming possible, which allow accurate treatment of both 
optically thick and thin pho ton fields in the "instant light" (nonrel- 
ativistic) approximation ([Haves & Normanl 120031 : iGonzalez et al.l 
l2007l:|jiang et al. 2012; Davis et al. 2012). However, as we discuss 
below, GR radiation transfer codes still continue to rely on flux- 
limited diffusion or the Eddington approximation. 

At super-Eddington accretion rates, L/Leaa > 1, accretion 



flows again become radiatively inefficient. Here, the optical depth 
is so large that the photon diffusion time from the disk inte- 
rior to the photosphere becomes longer than the accretion time. 
As a result, most of the photons are advected with the gas into 
the B H, leading to a "slim accretion disk" dAbramowicz et all 
1988). This important but poorly understood accretion state may 
be responsible for much of the SMBH mass growth in the Uni- 
verse; for instance, it might explain the paradox of having IO'Mq 
BHs alrea dy at z > 6, when the Universe was less than 1 
Gvr old dBarth et al I liooll; Iwillott et all 120051; iFan et al.l 120061; 
IWillott et al.l l20ld : iMortlock et al.l 1201 ll) . Super-Eddington ac- 
cretion may also apply to ultra-luminou s X-ray sources (e.g., 
IWatarai et al.l200ll : lKawashima et alj2012l) . Clearly, to understand 
the accretion physics in these systems, radiation MHD models 
that self-consistently couple gas, radiation and magnetic fields, are 
crucial. Such models ar e also important fo r measuring BH spins 
dMcClintock et alj|201 ll : IStraub et al.ll201 lh . calculating the radia- 
tive e fficiency of super-Eddington accretion disks (e.g., ISad owski 
2009), and understanding large-scale c osmological feedback by 
radiation-driven outflows from AGN (see lFabianl20 1 2| f or an obser- 
vational review). Some efforts have already been made on attack- 
ing these important problems, fn particular, using a non-relativistic 
code and flux-limited diffusion, super-Eddington a ccretion flows 
have been simulated and their spectra computed (Ohsugaetal 



1 However. Il5ibi et alj J20T2I) have shown that the accretion rate of Sgr A* 
is near the limit of the regime where radiative cooling may be important. 



| 2003|: IOhsugall2006i; lOhsuga & Mineshigdl2oTTI : [Kawashi ma et al. 
l2012h . 

Until a few years ago, all radiation hydrodynamic and MHD 
simulations of accretion disks were run with non-relativistic or 
special relativistic codes. However, physics around BHs must be 
studied in GR because of the the strong gravity involved. Re- 
cently, pro gress has been m ade in implementing radiation into 
GR codes. iFarris et all d2008l) developed a formalism for incorpo- 
rating radiation under the Eddington approximation in a conser- 
vative MH D code. This method has been implemented in other 
codes dZanotti etalj|201ll ; IFragile et alj|2012l) . but all these codes 
have problems because of the stiffness of radiative source terms 
at large optical depths. These stiff terms require implicit treat- 
ment, which is p rohibitively complicated in curved space-time. 
lRoedigetal1d2012h . extended the method by applying an implicit- 
explicit Runge-Kutta numerical scheme. However, the authors 
again relied on the Eddington approximation, which cannot handle 
optically thin flows ac curately. A more advanced method has been 
recently described by Shibata & Sekiguchi (2012]), who employ a 
truncated moment formalism, similar to our method, for neutrino 
transport in GR. 

In the present paper, we describe a simple approach for includ- 
ing radiation transport in GR codes. The method makes use of a 
key simplification that is intrinsic to GR applications. Normally, in 
non-relativistic codes, hydrodynamic or MHD signal speeds, which 
determine the time step via the Courant condition, are much slower 
than the speed of light. Correspondingly, the time step that one 
uses to evolve the fluid equations is much longer than the light- 
crossing time across a cell. This is a problem when one wants to 
simulate radiation hydrodynamical systems, especially in optically 
thin regions, where one is faced with a large mis-match between the 
characteristic speeds of the fluid and the radiation. Either one must 
limit the time step to the light-crossing time, which prohibitively 
increases the computational cost (as the fluid dynamics evolve on a 
much slower time scale), or one must handle all the radiation terms 
via an implicit method. The latter inevitably couples neighboring 
cells and makes the code very complicated. Moreover, it does not 
easily generalize to curved space-time. 
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In contrast, a GR code is applied only in relativistic space- 
times near BHs and neutron stars. The time step in simulations is 
generally set by applying the Courant condition to the smallest grid 
cell, located at the innermost radius, where the flow is relativistic. 
Thus, the normal time step in a pure GR hydro or GRMHD code is 
already limited by the speed of light. Therefore, including radiation 
as an extra relativistic fluid is fairly easy. In particular, the advec- 
tion terms in the radiation equations are no more difficult to com- 
pute than the corresponding terms in the fluid and magnetic field 
evolution equations. Nor does the time step need to be adjusted in 
any way. Because of this large simplification, the advective radia- 
tion operator can be treated via a standard explicit approach, just as 
one handles the corresponding hydro and MHD terms. 

Of course, the interactions between radiation and gas via emis- 
sion, absorption and scattering introduce their own time scales. 
These can sometimes be very short, requiring implicit handling. 
However, these interactions are local and can be handled via a lo- 
cal implicit scheme. This is a great simplification. It means that one 
can do implicit evolution independently in each grid cell, without 
coupling to other cells. Therefore, there are no space-time curvature 
effects to contend with as in other more sophisticated multi-cell im- 
plicit schemes in GR. 

In the work described here, we have implemented the above 
approach, closing th e radiation moment equations using the 
Ml closure scheme dLever more 1985 iDubroca &Feugeaslll999T; 
Gonzalez e t alj|2007h . Ml closure allows a limited treatment of 
anisotropic radiation fields and works well in both optically thick 
and thin regimes. We have implemented our method in a GR radi- 
ation hydrodynamics (GRRHD) code KORAL. The structure of the 
paper is as follows: In Section|2]we introduce the equations, in Sec- 
tion[3]we describe the numerical algorithm used by KORAL, in Sec- 
tion|4]we present a set of test problems which validate the scheme, 
and in Sectionf5]we discuss possible applications of the code. 



2 EQUATIONS 

2.1 Conservation laws 

A pure hydrodynamic flow is described by the following conserva- 
tion laws, 



(1) 

(2) 



where p is the gas density in the comoving fluid frame, vt is the 
gas four-velocity as measured in the "lab frame", and 7f is the 
hydrodynamical stress-energy tensor in this frame, 



7^ =(p + U + P)lFUy + POX, 



(3) 



with if and p representing the internal energy and pressure of the 
gas in the comoving frame. 

In the case of radiation hydrodynamics, it is conve- 
nient to introduce the r adiation stress-energy tensor 7^ (e.g., 
Mihalas & Mihalai] 1 19841) . and to replace the second equation 
above with the more general conservation law, 



(2? + *S)».=0. 



(4) 



The radiation stress-energy tensor in an orthonormal frame com- 
prises various moments of the specific intensity 7 V , e.g., in the fluid 
frame it takes the following form, 



R 



E F' 



(5) 



where the fluid-frame quantitie^] 

E = f T v dvdQ, 



F' 



f 
f 



I v dvd£lN, 
IdvdflNN 



(6) 
(V) 
(8) 



are the radiation energy density, the radiation flux and the radiation 
pressure tensor, respectively, and N' is a unit vector in direction x' . 

The fluid frame radiation stress tensor R is related to the tensor 
R defined in the locally flat non-rotating frame, o r the zero-angular 
momentum frame (ZAMO) l lBardeen et al.lll972h . by 



where A is the Lorentz boost, 

y 



A(«) : 



yv J 6' J + 



yv 

v l v J (y—Y) 



(9) 



(10) 



7 = u', v' = u'/u', and if is the four-velocity of the gas as measured 
by the locally non-rotating observer. 

Quantities in the ZAMO frame (denoted with tildes) are re- 
lated to those in the lab frame by tensors created from the compo- 
n ents of the correspond ing tetrads of the ZAMO, e£ and ej, defined 
in lBardeenetaT] ( ll972h . The radiation tensor transforms as 



R" y 
RK" 



while four- vectors transform as 

vf = e*iu". 



(11) 
(12) 

(13) 
(14) 



The conservation law © may be rewritten with the help of the 
radiation four-force density G" as 



(T{%, = Gy, 

(R" v ),, = ~Gy, 

where G" is given by <Mihalas&Mihak3 [T984h. 

G v = f t*r v / v -7z v )dvdQN\ 
which takes a particularly simple form in the fluid frame. 



| K{E-AnB) 



(15) 



(16) 



(17) 



Here, B = <tT 4 /jt is the integrated Planck function corresponding 
to the gas temperature T, a is the Stefan-Boltzmann constant, Xv 
and r)y denote the frequency-dependent opacity and emissivity co- 
efficients, respectively, while k and x are the frequency integrated 
absorption and total opacity coefficients, respectively. In Section|4] 
we occasionally refer to the scattering opacity k cs , which is related 
to k axidx by 



X = K + K es . 



(18) 



The four-force C may be transformed between frames as described 
above, e.g., 



(19) 



2 Throughout the paper, "widehats" denote quantities in the fluid frame and 
"tildes" denote quantites in the ZAMO frame. 
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The rest mass conservation equation (Q] and the energy- 
momentum conservation equations d!5t may be w ritten in a co- 
ordin ate basis in the following conservative form dGammie et al.l 
120031) . 

d,(^pu') + d,(^pu') = 0, (20) 

a,( v=in) + v=in) = + ^Fsg v , (2D 

d^^R'^ + dii^K) = V=i^- V=gG v . (22) 

This formulation has a drawback for numerical computations: the 
terms involving Christoffel symbols on the right, when calculated at 
cell centers, will not balance the corresponding spatial derivatives 
on the left (approximated under a given reconstruction scheme). 
This is true even for particularly simple situations such as constant 
gas or radiation pressure, and can lead to catastrophic secular er- 
rors. To solve this issue we can either modify the values of the 
Christoffel symbols suitably (Appendix A of McKinney et al. 2012) 
or we can reformulate the equations so as to avoid the problem. We 
choose the second approach and make use of the following equa- 
tions, 



written in a compact form as 



d,(pu') + di(pu') 



pu 



3,(V=gX 



(23) 
(24) 



d,(T' v ) + d t (Tl) = r t ri - -^d,( V=g) + G„, 

d,(R' v ) + = R^Tl - -^La,( V=D - G r , (25) 

where we assumed that the metric is static and moved its determi- 
nant out of the derivatives on the left. In this formulation, the two 
terms that are expected to cancel each other both appear as source 
terms and therefore are calculated at the same locatior0. 



R" v = -Eitiu R +-Ef\ 

3 3 s 



(27) 



where = (1,0,0,0) and gf" is the contravariant metric tensor, 
which in this frame is given by the flat space Minkowski metric. 
Since equation \21\ is in a covariant form, it is also valid in the 
lab frame (as well as all other frames), with if R being the four- 
velocity of the radiation rest frame as measured by an observer in 
the lab frame. Note that, regardless of which frame one works in, 
the quantity E should be interpreted as the radiation energy density 
as measured in the radiation rest frame. 

Each time step in the numerical integration in any particular 
cell gives an update to the "time" row of the radiation tensor, R' v , 
for that cell in the lab frame. Thus, we obtain numerical values of 
these four particular components of the tensor. According to equa- 
tion d27b , the full tensor R 1 " is a function of 5 numbers, E, if R , 
though only four of these are independent since the norm of the 
four- velocity if R is equal to -1. Hence, we can use the four given 
tensor elements to solve for the four unknowns and thereby com- 
pute the full matrix. Below we give an algorithm for doing this 
analytically. 

Consider the quantity R m R' v which can be expressed as (equa- 
tionl27t. 

RVRlv = 1 £ 2 [i6(4)X^ + K«£g<* + Au< R ulg'» + g'V] • 

(28) 

If we contract this with and use the following results, 

= -1, (29) 
g/ivKs™ = Uv,rS"' = u 'r> (30) 
S' v g' v = g", (31) 



ft V 
8fiV U R U R 



we obtain 



2.2 Closure scheme 

To close the above set of equations we need a prescription to com- 
pute the second moments of the angular radiation intensity distri- 
bution. Specifically, we need a prescription to write down the full 
radiation stress tensor R^ v knowing only the radiative energy den- 
sity R" and the fluxes R". 

The simplest approach, which corresponds to assuming a 
nearly isotropic radiation field, is the Eddington approximation, 
which in the fluid frame gives 



pi = Le6 u . 



(26) 



However, the assumption of isotropic specific intensity is good only 
in optically thick media. In many astrophysical applications we are 
interested in radiation that escapes from the photosphere to infinity, 
for which we n eed a better c losure scheme. 

Following iLevermorel J 19841) . we assume that the radiation 
tensor is isotropic and satisfies the Eddington closure, not in the 
fluid frame, but in the orthonormal "rest frame" of the radiation. 
The latter is defined as the frame in which the radiative flux van- 
ishes. Thus, in this frame, we assume that R" = E, R" = E/3, and 
all other components of R are zero. This leads to the Ml closure 
scheme. 

In the radiation rest frame, the radiation stress tensor can be 



g, v W = -^ 2 K) 2 + i£V. 



(32) 



The left-hand side is computable from the four given tensor ele- 
ments and the right-hand side involves two of the unknowns: E, u' R . 
We also have the following expression for R", 



R" = -E(u' R ) 2 



+ 3 V. 



(33) 



which again involves the same two unknowns. Thus, we can solve 
equations d32t and {33} to obtain E and u' R (it reduces to a quadratic 
equation). It is then straightforward to calculate the remaining u' R 
from the other time components of equation ( 127) and to calculate 
the entire radiation stress tensor. 

For flat space time, the above f ormulation reduces to the 
standard formulae (Levermorel 1 19841 : iDubroca^ fe Feugeas; 1 19991: 
iGonzalez et alj|2007f) . For instance, the radiation pressure tensor 
P'i in the fluid frame has the form, 



2 2 \fp) 



(34) 



where /' = F' IE is the reduced ra diative flux and £ is the Eddington 
factor given by dLevermorell 1 984b . 



3+4/7; 



5 + 2V4-3/'/- 



(35) 



3 In principle, it is sufficient to move the determinant out of the spatial 
derivatives in r- and S- components of equations 12 It and (22) . 



In the extreme "optically thick limit", F' ~ 0, and we find 
/' = 0, f'fi = and <f = 1/3, which corresponds to the correct 
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answer, viz., the Eddington approximation, 



1/3 






1/3 
1/3 



E. 



(36) 



In the opposite extreme "optically thin limit", F l = E, i.e., a uni- 
directional radiation field directed along the x-axis, we have /' = 
S\ > f'fi = 1 an d f = 1/3, which gives 



ptj 



1 






(37) 



This corresponds to an intensity distribution in the form of a Dirac 
^-function parallel to the flux vector, which is again the correct an- 
swer. The Ml closure scheme thus handles both optical depth ex- 
tremes well and it is found to be fairly good at intermediate optical 
depths as well. 

As explained, the Ml closure scheme assumes that radiation is 
isotropic in the radiation "rest frame". The stress tensor in an arbi- 
trary frame is obtained by applying a Lorentz boost to the isotropic 
"rest frame" tensor. As a result, only one direction, the direction 
of the boost, is distinguished. In other words, the specific inten- 
sity is always symmetric with respect to the mean flux. The Ml 
closure scheme is thus expected to be only approximate when mul- 
tiple sources of light are involved (see Section 14.5) . In problems 
involving accretion disks, which are the primary area of interest of 
the present authors, highly anisotropic configurations with multiple 
beams are not very common, and the Ml scheme is probably ade- 
quate. In any case, Ml closure will provide a significantly superior 
treatment of radiation in the optically thin regions near and above 
the disk photosphere, compared to the Eddington approximation or 
flux-limited diffusion. 



3 THE KORAL CODE 

The scheme described in this paper has been implemented into a 
GRRHD code KORAL which solves equations J23M25b in an arbi- 
trary metric. The code uses a finite difference scheme with eit her 
linear slope-limited reconstruction (Kurganov & Tadmor 2000l) or 
fifth-order non-linear monotonizing filter reconstruction (MP5, see 
ISuresh & Huvnhlll999l ; iDel Zanna et alj|2007h . The fluxes at the 
cell faces are calculated using the Lax-Friedrichs scheme. The 
source terms are applied at the cell centers and the time stepping 
is performed u sing the optimal Runge-Kutta method of third order 
(Shu & Oshedll984l) . The vector of conserved quantities is (Sec- 
tionl3~5l 



U = \pu\T\ + pu' ,T\,R' t ,R\l, 
while the primitive quantities are, 

P= \p,u,v?lvt,E>F\ 



(38) 



(39) 



Conversion from conserved to primitive quantities is described in 
Section [3~4l while the algorithm itself is described in the next one. 



3.1 The algorithm 

During each sub-step of the Runge-Kutta time integration, the code 
carries out the following steps in the given order: 

(i) The vector of conserved quantities in each cell is used to cal- 
culate the primitive quantities at the cell center (see Section [J!4t . 



(ii) Ghost cells at the boundaries of the computational domain 
are assigned primitives appropriate to the boundary conditions of 
the particular problem of interest. 

(iii) For each cell, the maximal characteristic left- and right- 
going wave speeds are calculated, following the algorithm de- 
scribed in Section [37H 

(iv) For each dimension, primitives are interpolated using the 
chosen reconstruction scheme (linear slope-limiter or non-linear 
monotonizing filter) to obtain their left- and right-biased values at 
cell faces: P L and P R . 

(v) From P L and P R , left- and right-biased fluxes T L and T R are 
calculated at cell faces. 

(vi) The flux at a given cell face is calculated using the Lax- 
Friedrichs formula, 



T = l -{T R + T L -a(U R -U L )), 



(40) 



where a is the maximal absolute value of the characteristic speeds 
at the centers of the two cells on either side of the face, and U R and 
Ul are the conserved quantities calculated at the cell face based on 
Pr and P L . 

(vii) The advective time derivative is calculated using an unsplit 
scheme, 



dU_ 

dt (adv) 



n-n 

dx 1 



dx 2 



n-n 

dx 3 ' 



(41) 



where dx 1 denotes cell size in the direction i. Note that all prim- 
itives, including the radiation density E and radiation flux F*, are 
treated identically as far as the advective term is concerned. 

(viii) The geometrical source terms, viz., all terms on the right 
hand sides of equations l !23M25t except the radiation four-force 
density ±G„, are calculated at cell centers to give the corresponding 
time derivative dU/dt^ eo) . 

(ix) The advective and geometrical operators are used to update 
the conserved quantities according to 



AU : 



dU dU 
— + — At 

at (adv) at (gco)/ 



(42) 



That is, all these terms are treated in an explicit fashion. 

(x) The updated vectors of conserved quantities are used to cal- 
culate the corresponding updated primitive quantities at cell cen- 
ters. 

(xi) Finally, the remaining terms, viz., those involving the four- 
force density C, are handled implicitly using the method described 
in Section [33] This results in a final update of the vector of con- 
served quantities at each cell center. 

3.2 Characteristic wavespeeds 

The Lax-Friedrichs scheme requires knowledge of the maximal 
characteristic wave speeds of the system (a in equation I40t. The 
hydrodynamical and the radiative components of equations (120b- 
d22b are coupled only through the radiative source term ±G V . There- 
fore, for the purpose of calculating the fluxes at cell faces and the 
advective time derivative, we are allowed to separate the hydro- 
dynamical and radiative wave speeds. We can calculate each sep- 
arately, and use it for evaluating its corresponding flux. Such an 
approach avoids excessive artificial numerical viscosity which ap- 
pears when the characteristic wavespeeds are not separated. The 
radiative wavespeed, as described below, never drops below cj V3. 
If such a high value was used in equation l |40| > for the hydrodynam- 
ical subsystem, it would result in strong artificial diffusion of the 
gas. 
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The hydrodynamical characteristic velocity in the fluid frame 
is the local speed of sound, 



Tp 



p + u + p 



(43) 



To get the left- and right-going wave speeds in the lab frame we 
transf orm this velocity in the standard way (e.g., iGammie et al.l 
l2003h . 

The evolution of the radiation field is described in the fluid 
frame by the following set of equations, 



d,E + diF' 
d,F' + dp' 



-G', 
-G J , 



(44) 



which is the non-relativistic limit of equations (122b in flat space- 
time. Given a closure scheme, it is possible to calculate the Jacobi 
matrices, 



dF/dE 
aP ij /dE 



dFjdF> 
dP'IdF' 



(45) 



whose eigenvalues give the wave speeds of interest in a given di- 
rection i. The maximal left- and right- going speeds can then be 
transformed to the lab frame following the same method as for fluid 
velocities. 

For the Ml -closure scheme, the equation for the eigenvalues 
is quartic and may be solved efficiently using standard algorithms. 
Another approach, which will likely improve code performance, is 
to precalculate tables of radiative wave speeds as a function of the 
reduced flux vector components F'/E (Gonzal ez et alj|2~007h . At 
present, KORAL uses the analytical approach. 

In the limit of large optical depths, the radiative energy den- 
sity, when decoupled from gas (e.g., for «« 1 but^ » 1), has a 
diffusion coefficient D given by (see Section l4~4l 



1 

D = — . 

3X 



(46) 



In this limit the distribution of radiative energy density should 
remain stationary (d/dt — > 0). On the other hand, the maxi- 
mal eigenvalue of the Jacobi matrices J' (equation 1451 ) is ± 1 / V3 
dGonzalez et al.|[2007h . If such large wave speeds are incorporated 
into a numerical scheme they will result in large, unphysical, nu- 
merical diffusion. To limit this effect, we modify the radiative wave 
speeds in the fluid frame according to 



3t' 



(47) 



3t> 



where a' R and a' L are the maximal right- and left-going radiative 
wave speeds in the fluid frame in the direction i, and r' = ^dx 1 is 
the total optical depth of a given cell in that direction. 

The smaller the characteristic wave speed in equation (14 It . the 
weaker the numerical diffusion. Thus, one may be tempted to set 
the wave speed to zero. However, the numerical scheme will then 
no longer satisfy the total variation diminishing (TVD) condition 
and the algorithm will be unstable. Our choice of the wave speed 
limiter (equation I47t is motivated by the fact that, for a diffusion 
equation of the form y j( = Dy xx , the maximum allowed time step 
for an explicit numerical solver is 



At ■ 



(.Ax) 2 
AD 



(48) 



This expression, combined with equation d46t . gives 
Ax 4 4 

— = = — , (49) 

At 3xAx 3t y ' 

which is the limiter introduced in equation ( 147) . 

3.3 Implicit treatment of radiative source terms 

It is well-known that, under some circumstances, e.g. for large op- 
tical depths, the radiative source terms ±G V in equations d21t and 
( 122b bec ome stiff, making e xplicit integration practically impossi- 
ble (e.g.. lZanottietalJl201l|) . We then need to treat these terms via 
implicit time integration. In principle, we could try to solve the 
whole system of partial differential equations d23)-d25) implicitly, 
as done for instance in non-relati vistic or special relativistic radia- 
tion h ydrodynamics codes (e.g., Kru mholz et alj ;2007; Jian g et all 
l2012h . However, this approach is very difficult in GR, where the 
curvature of spacetime makes the problem highly complicated and 
it is non-trivial to ensure that an implicit code is conservative. 

As already explained, our approach is to split the advective 
derivative operator from the radiative source terms operator. The 
former is applied explicitly in the usual way while the latter is han- 
dled implicitly. This approach is possible because the time step is 
already limited by the speed of light just from the fluid dynamics, 
so radiation advection is also guaranteed to be stable in an explicit 
scheme0The main advantage is that the radiation source terms ±G V 
are local, so they can be treated semi-implicitly and point-wise in 
the fluid frame, without having to deal with curvature effects. In our 
experience, this approach is both simple and robust. 

The radiative source term operator describes the action of the 
radiative four-fource G* on the energy and momentum density of 
the gas and the radiation. The corresponding equations are 



d,(T' v ) = G y , 

d,(R[.) = -Gy. 



(50) 
(51) 



In an explicit scheme, updates would be calculated very simply as 



of _ pt _ 

n v,(»+l) v,(n) — 



At G Vj( „), 
-At G vM , 



(52) 
(53) 



where the subscripts (») and (n + 1) denote values at the beginning 
and end of a time step of length At, respectively. This approach, 
though simple, is numerically very unstable whenever the terms in 
the force vector G v are large. Our scheme avoids the instability by 
computing the updates implicitly via 



R'„ 



V,(J1+1) 



. T< - 

' K(n) = 



At G v> („ + i), 

-At Gy_(„ + 1), 



(54) 
(55) 



i.e., using quantities at time (n + 1) rather than (;;) to compute the 
force vector on the right-hand side. It is well-known that this simple 
change has a profound effect on stability. 

KORAL solves equations d54t and \55\ numerically. Because of 
the symmetry of the problem, specifically, the right-hand sides of 
eqs. d54| l and d55t differing only by sign, the system of equations 
may be reduced to four non-linear equations, e.g., for R', n+V) - We 
use the standard Newton method to solve these equations and esti- 
mate the Jacobian matrix numerically. During each iteration we use 

4 Note that if we were to apply the code to non-relativistic problems, or 
situations involving weak gravity, there would be a large disparity in the 
advective time scales of fluid and radiation, and the code would no longer 
be efficient. 
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the current guess of R' . +1 . to calculate the corresponding values of 
V 

We then convert the current vector of conserved quantities to prim- 
itives (see Section [3T4l and calculate the radiative four force in the 
fluid frame, G y (n+l) . This is then boosted to the laboratory frame to 
obtain G v ,(„+i). 

The method described above is numerical and can occasion- 
ally fail. In Appendix lAl we describe an alternate, fully analytical, 
but approximate, method for applying the radiative source opera- 
tor. KORAL uses that algorithm as a failsafe backup whenever the 
numerical scheme descibed in this section fails. 

3.4 Conversion between conserved and primitive quantities 

In conservative GR numerical codes it is necessary to convert con- 
served quantities to primitives at least once per sub-time step (twice 
in the case of our algorithm). In our problem, the hydrodynamical 
and radiative variables decouple, so the conversion may be done 
separately for each. 

The conversion of hydrodynamical quantities is performed i n 
the usual manner (e.g.. iNoble et alj|2006t iDel Zanna et al.ll2007l) . 
We take the explicit forms of the conserved quantities (equation^, 
use the equation of state p = (T - l)u, the four- velocity normal- 
ization u^u^ = -1, and combine all the terms into a non-linear 
equation for the internal energy density u. We solve this equation 
numerically by the Newton method using the value of u from the 
end of the previous time step as the initial guess. 

The radiative conserved variables (R' ) may be easily con- 
verted to the radiative primitives (E, F '). First, one has to calculate 
all the components of the radiation stress tensor in the lab frame 
Rf v following the algorithm described in Sect. 12.21 This tensor is 
then boosted to the fluid frame (Sect. |2~Tt , 

R» v = AgH&AJHOajSj**. (57) 

The fluid-frame radiative energy density E and fluxes F' are then 
given by, 

E = R", (58) 
F' = R". (59) 

3.5 Implementation notes 

(i) The mass conservation equation ( |23t and the gas internal en- 
ergy conservation l aw, i.e., the t compo nent of equation {24}, are 
aggregated to give dGammie et al .2003) 

d, {t; + P u<) + d, (t; + P iA = z*r* - Tv + J^' a,( y=J) + g„ (60) 

which replaces the t component of equation d24t . Then, (T' t + pu') 
becomes the relevant conserved quantity, which reduces in the non- 
relativistic limit to T\ + pu' — > -it. 

(ii) In cold relativistic flows, where u <k p, the numerical accu- 
racy is not suficient to evolve the internal energy reliably. As a re- 
sult, negative internal energy densities may be occasionally found. 
Whenever this happens, we recalculate the gas properties by evolv- 
ing the gas entropy per unit volume, 

^(rV 0S # (61) 



Table 1. Model parameters for relativistic hydrodynamical shock tube test 





Left state: 


Right state: 


r 


p p u x 


p p if* 


5/3 


10.0 13.33 0.0 


1.0 10~ 8 0.0 



That is, we compute the entropy at the end of the last successful 
time step and evolve it according to (5 uf) ;ll = 0. 



4 TEST PROBLEMS 
4.1 Relativistic shock tube 

To validate the hydrodynamics part of our code, w e tested it us- 
ing the relativistic shock tube problem introduced bv lHawlev et al.l 
( 1984). The modeled system consists of two states (left and right), 
separated initially by a membrane. The gas on the left is dense and 
hot, while that on the right is rarefied and cool (Table [TJ. At time 
t = the membrane is removed and the hot gas pushes the cool 
gas to the right causing a shock to travel to the right. Meanwhile, 
a rarefaction wave moves at the local sound speed back to the left. 
An analytical solution for this shock tube problem may be obtained 
by solving the appropriate jump conditions. 

In Figure [T] we show our numerical solutions at time t = 50, 
which may be compared with the exact so lution (black dashed line) 
obtain ed using an exact Riemann solver dGiacomazzo & Rezzollal 
2006). The upper panel shows profiles of density obtained using 
three different reconstruction schemes: linear slope limiter with 
6=1, linear slope limiter with 6 = 20 and fifth-order monotoniz- 
ing filter, MP5. The rarefaction wave and the plateau are resolved 
equally well in all three schemes. However, there are differences in 
the post shock region. The linear 9=1 scheme (red lines) is most 
diffusive and smooths the edges of the postshock region, while the 
MP5 scheme (blue) somewhat underestimates the density here. The 
lower panel shows the velocity. This is reproduced well in all three 
schemes, though the MP5 scheme produces low-level oscillations 
near the edge of the rarefaction wave. 

This test shows that all the reconstruction schemes currently 
implemented in KORAL work reasonably well on relativistic hydro- 
dynamic shocks. The MP5 scheme is most accurate in smooth re- 
gions, but it is also prone in some circumstances to give unphysical 
oscillations. The 9=1 linear scheme is most diffusive, and at the 
same time is also the most stable. In the following tests, if not stated 
otherwise, we use the linear reconstruction scheme with 9 = 1.5. 



4.2 Polish doughnuts 

To test the ability of the code to handle multi-dimensional hydrody- 
namics in curved s pacetime, we set up analy tical equilibrium torii 
(Polish doughnuts, lAbramowicz et al.ll 1 9780 in the Schwarzschild 
metric as initial conditions and let the code evolve to a numerical 
equilibrium configuration. For the analytical model, we assume a 
constant specific angular momentum, I = —u^ju t = constant. From 

5 The 9 parameter in the generalized minmod lim iter determines the diffu- 
sivity of the scheme (Kurganov & Tadmor 2000); 8 = 2 corresponds to the 
MC (monotonized central) scheme and 6=1 corresponds to the MINMOD 
scheme. 
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Table 2. Model parameters for Polish doughnuts tests 



Figure 1. The relativistic shock tube test solved on 500 grid points using 
three reconstruction algorithms: slope-limited linear with 6=1 (most dif- 
fusive, red), 8 = 2 (green), and MP5 (blue). The black dotted line shows the 
exact solution for comparison. 



the condition ifu p = -1, it follows that 

u; 2 = -g" + 2ig" f ' - fg* 



(62) 



We choose the specific internal energy at the inner edge of the torus, 
m, in , which determines the radius of the inner edge of the torus, and 
we th en calculate the fluid enthalpy, h = o+u+v (e.g. jHawlev et ail 
11984) , 



h : 



(63) 



Using an equation of state p = Kp T (where the constant k determines 
the entropy of the torus gas), we obtain 



(h-l) (T- 1) 



i/Cr-i) 



(64) 



(65) 



= 0, v* 



uf/u', and choose 



We set the initial velocity to v r 

r = 4/3. 

Table[2]gives parameter values corresponding to three models 
that we ran to test the code. Models 1 and 2 have the same value of 
the specific angular momentum, £ = 4.5. Model 1 has u,- m = -1, 
corresponding to a torus inner radius r m = 8, while Model 2 has 
t'r.in = -0.98, corresponding to r m = 10. The specific angular mo- 
mentum for the third model, £ = 3.77, lies in between the Keplerian 
values of £ at the marginally stable and marginally bound orbits 
(t ms < £ < l„,b)- Therefore, this torus has a cusp (self-crossing 
equipotential surface) near the marginally bound orbit (r m t = 4). 
All the equipotential surfaces outside the body of the torus (defined 
by the critical surface producing the cusp) are open and continue 
into the BH. Therefore, the assumption of zero poloidal velocity 
cannot be applied to this region. 

The above three torus models were simulated on a 50x50 grid 
in Boyer-Lindquist {r-0) coordinates, linearly spanning the range 
r = 3 - 27.8 and 9 = 0- n/2. At the spin axis and the equatorial 
plane, reflection boundary conditions were assumed. The analytical 
solution was imposed as boundary conditions in the ghost cells out- 
side r = 27.8. At the inner edge of the grid r = 3, outflow boundary 



Model 


I 


"t.in 


K 


1 


4.5 


-1.0 


0.03 


2 


4.5 


-0.98 


0.03 


3 


3.77 


-1.0 


0.03 



conditions were implemented. Linear reconstruction with G = 1.5 
was used. 

Figure [2] compares the relaxed, stationary solutions obtained 
by running the code with the corresponding analytical Polish 
doughnut solution. Colors and orange solid lines show the numeri- 
cal solution (density), while dashed green lines show the analytical 
model. Vectors indicate the poloidal component of velocity. 

The top panel in Figure ^corresponds to Model 1, where the 
analytical solution extends from r = 8 to the outer boundary. The 
analytical and numerical contours agree very well. The only visible 
discrepancy is in the low-density region near the torus inner edge, 
where there is some numerical dissipation. 

The middle panel in Figure [2] corresponds to Model 2, where 
the analytical torus is entirely confined within the grid: r m = 10, 
r out = 27. Again, the numerical solution agrees very well with 
the analytical solution except for the innermost region of the torus 
where the numerical solution is slightly stretched inward. 

The bottom panel in Figure [2] shows Model 3. The analytical 
contours plotted correspond to open equipressure surfaces which 
intersect the BH horizon. Non-zero poloidal velocities are thus ex- 
pected in this model and the analytical solution is not fully consis- 
tent. Nevertheless, there is a close similarity between the numerical 
and analytical solutions, which means that the code is able to rep- 
resent even this torus quite well. 



4.3 Radiative shock tubes 

The previous two tests did not involve radiation. For our first test 
with radiation, we s et up a number of radia tive shock tube prob - 
lems as described in iFarris et all d2008l) and iRoedig et ail d2012h . 
These shock tubes are similar to the pure hydrodynamic shock tube 
described in Section l4"Tl i.e., the system begins with gas in two dif- 
ferent states (left and right), separated by a membrane. The mem- 
brane is removed at t = and the system is allowed to evolve. 
The difference here is that the evolution is described by the full set 
of equations fl23|l-(|25t. In addition, the left- and right-states of all 
the tests except test No. 5 are set up in such a way th at the shock 
asymptotically becomes stationary (see Appendix C of lFarris et al.l 
120081) . 

Table [3] lists the parameters describing the initial states of 
seven test problems that we have simulated. The scattering opac- 
ity in all the tests is set to zero, so x = k (equation 1 18t. The value 
of the radiative constant <x in code units is given in the table. All 
the tests were solved on a grid of 800 uniformly spaced points. For 
consistency with previous work by other authors, the Eddington 
approximation was used (equation|26). 

Figure [3] shows the numerical (solid) and analytical (dashed 
lines) solutions for radiative shock tube problem No. 1, which cor- 
responds to a non-relativistic strong shock. The panels show (top to 
bottom) density, proper velocity, fluid-frame radiative energy den- 
sity and fluid-frame flux, and may be directly compared (but for 
the flux which the other authors plot in the lab frame) to the cor- 
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Table 3. Radiative shock tubes 



Test Left state: Right state: 



No. 


r 


(T 




k/p 


P 


p 




£ 


P 


P 


u x 


E 


1 


5/3 


3.085 


10 9 


0.4 


1.0 


3.0 X IO" 5 


0.015 


1.0 x 10~ 8 


2.4 


1.61 x 10~ 4 


6.25 x IO" 3 


2.51 x IO" 7 


2 


5/3 


1.953 


10 4 


0.2 


1.0 


4.0 X IO" 3 


0.25 


2.0 x 10~ 5 


3.11 


4.512 x 10~ 2 


8.04 x 10~ 2 


3.46 x 10~ 3 


3a 


2 


3.858 ■ 


io- 8 


0.3 


1.0 


6.0 x IO 1 


10.0 


2.0 


8.0 


2.34 x IO 3 


1.25 


1.14 x IO 3 


3b 


2 


3.858 ■ 


10- 8 


25.0 


1.0 


6.0 x IO 1 


10.0 


2.0 


8.0 


2.34 x IO 3 


1.25 


1.14 x IO 3 


4a 


5/3 


3.470 


IO 7 


0.08 


1.0 


6.0 x 10~ 3 


0.69 


0.18 


3.65 


3.59 x 10~ 2 


0.189 


1.3 


4b 


5/3 


3.470 


10 7 


0.7 


1.0 


6.0 x IO -3 


0.69 


0.18 


3.65 


3.59 x 10~ 2 


0.189 


1.3 


5 


2 


3.858- 


10~ 8 


1000.0 


1.0 


6.0 x IO 1 


1.25 


2.0 


1.0 


6.0 x IO 1 


1.10 


2.0 




Figure 2. Equilibrium torii (Polish doughnuts) around a Schwarzschild BH 
after numerical evolution up to time t = 1000M. Colors and arrows show 
density and poloidal velocity, respectively. The green dashed lines show 
iso-density contours of the analytical solutions. Contours for the numerical 
solutions are shown with orange solid lines. The parameters of the three 
model torii are given in Table [2] 



responding figures and analytical solu tions of dFarris et al.ll2008l : 
IZanotti et al]201ll ; lFragile et alj|2012l) . The agreement is good. 

Figure [4] shows results for radiative shock tube test No. 2, 
which corresponds to a mildly relativistic strong shock . Again, the 
agree ment between the numerical and semi-analytical dFarris et al.l 
2008) profiles is good, except for a slight smoothing of the numer- 
ical profiles at the position of the shock (see the bottom panel). 
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Figure 3. Results obtained for radiative shock tube test No. 1. From top to 
bottom, the panels show the profiles of gas density, proper velocity, fluid- 
frame radiative energy density and fluid-frame radiative flux. The numeri- 
cal solution obtained with KORAL is indicated by solid lines and the semi- 
analytical solution given in Farris et al. (2008) by dotted lines. 
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Figure 4. Same as Figure|5]but for radiative shock tube test No. 2. 
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Figure 5. Similar to Figure [3]but showing results for radiative shock tube 
tests No. 3a (solid blue: numerical solution, dotted green: exact semi- 
analytical solution) and No. 3b (solid red: numerical, dotted magenta: semi- 
analytical). 



4 F 




-5 5 10 

X 



Figure 6. Similar to Figure [3] but showing results for radiative shock 
tube tests No. 4a (solid blue: numerical solution, dotted green: semi- 
analytical solution) and No. 4b (solid red: numerical, dotted magenta: semi- 
analytical). 



Figure [5] shows results corresponding to radiative shock tube 
tests No. 3a and 3b. These are strongly relativistic shocks with up- 
stream u x = 10. Test No. 3a corresponds to shock tube test 3 of 
Farris et al while test 3b is the optically thick version of th e 

same test which was proposed and solved by Roedig et alj d2012h . 
These two tests verify that the code is able to resolve a highly rel- 
ativistic wave in two very different optical depth limits. In both 
cases, the numerical solution reaches a steady state and closely 
follows the corresponding semi-analytical solution. The only no- 
ticeable difference is that in the high-opacity case the discontinuity 
in the numerical solution is less steep than in the exact analytical 
profile. This discrepancy measures the amount of numerical dissi- 
pation introduced by the algorithm. 

Figure|6]shows results for radiative shock tube tests No. 4a and 
4b. These tests correspond to radiation pressure dominated mildly 
relativistic waves. Tes t 4b is the op t ically thick version of test 4a 
that was proposed by iRoedig et al.l d2012h . In both tests, the nu- 
merical solution reaches a stationary state and agrees well with the 
semi-analytical solution. Note that the values of the opacity coeffi- 
cient k in t ests 3b and 4b are t he maximum values that the numerical 
scheme of Roed ig et alj J2012h could handle. The algorithm imple- 
mented in KORAL has no such a limitation. We could increase k to 
much larger values and the scheme would remain stable. 

Figure[7]corresponds to radiative shock tube test No. 5. This is 
the only test that does not asym ptote to a stationary solution. This 
test was proposed and solved bv lRoedig et alj ( 120121) and represents 
an optically thick flow with mildly relativistic velocities. The left- 
and right-states are identical except that they have different veloci- 
ties. As a result, two shock waves propagate in opposite directions. 
This test does not have an analytical solution. Howev er, by com- 
paring our numerical solution with that presented in Roe dig et alj 
we confirm that our scheme performs satisfactorily. 




-20 -15 -10 -5 5 10 15 20 
X 



Figure 7. Same as Figure[3]but for radiative shock tube test No. 5. There is 
no analytical solution available for this problem. 

4.4 Radiative pulse 

We now test the ability of our scheme to handle the evolution of a 
radiation pulse in the optically thick and optically thin limits. We 
start with the optically thin case. We set up a Gaussian distribution 
of radiative energy density at the center of a 3D Cartesian coordi- 
nate system. The pulse radiative temperature is set according to, 

r rad = (^)" 4 = T (l + lOO^ 2 ^' 2 ) , (66) 

with T = 10 6 , w = 5.0 and the radiative constant a = 1.56x 10~ 64 . 
We assume zero absorption opacity (k = 0) and non-zero but negli- 
gible scattering opacity (/c cs = 10~ 6 ). The background fluid field has 
constant density p = 1 and temperature T = Tq. We solve the prob- 
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lem in three dimensions on a coarse Cartesian grid of 51x51x51 
cells using the linear reconstruction with 8 = 2. 

The initial pulse in radiative energy density is expected to 
spread isotropically with the speed of light (optically thin medium) 
and to decrease inversely proportionally to the square of radius (en- 
ergy conservation). Such behavior is visible in Fig. [8] showing the 
radiative energy distribution in the z = plane (top panels) and 
its cross-section along y = z = (bottom panels). The orange cir- 
cles in the top set of panels show the expected size of the pulse. It 
is clear that the propagation speed of the pulse is consistent. This 
problem was solved on a relatively coarse Cartesian grid. This re- 
sults in deviations from the perfectly spherical shape — the radia- 
tive energy density of the pulse along the axes is higher than along 
the diagonals. This effect is reduced by choosing larger resolution 
or a more suitable grid (e.g., spherical). The bottom set of panels 
in Fig. [8] shows the profiles of the energy density along the x-axis. 
The black dotted lines show the expected rate of energy decrease 
with increasing distance from the center. The numerical solution 
perfectly follows this trend. 

To test the optically thick limit we choose to set up a similar 
pulse but this time planar instead of a point-like, i.e., according to, 

Tmd = ( — J = T (l + 100e-*l*) . (67) 

This time we set the scattering opacity to /c cs = 10 3 and solve 
the problem as one-dimensional on 101 grid points distributed uni- 
formly between x = -50 and x = 50 with periodic boundary con- 
ditions in y and z. The total optical depths per cell and per pulse are 
therefore t = 10 3 and r = 10 4 , respectively. 

In the optically thick limit the evolution of such a system is 
described by a diffusion equation, 

d,E = —d m E, (68) 
3* 

which can be derived from the non-relativistic limit of Eq. [25] as- 
suming the time derivative of the x component of the flux vanishes. 
An initially Gaussian pulse of radiative internal energy will there- 
fore diffuse according to the value of the diffusion coefficient 1/3^. 

In Fig.[9]we plot profiles of the radiative energy at various mo- 
ments (solid lines) and compare them to the exact solution of Eq.[68l 
(dotted lines). The numerical solution diffuses slightly faster due to 
the additional numerical dissipation introduced by the schem^|. At 
later times this difference becomes insignificant. 

4.5 Shadow test 

Here we test the ability of the Ml closure scheme, as incoiporated 
in KORAL, to resolve shadows. We set up a blob of dense, optically 
thick gas in flat spacetime, surrounded by an optically thin medium, 
and we illuminate this system. 

We start with a single source of light imposed on the left 
boundary. We solve the problem in two dimensions on a 100x50 
grid, with the density distribution set to be 

P = Po + (Pi - Po) e V^W , (69) 
where p = 10~ 4 , p\ = 10 3 and w = 0.22. The gas temperature is 



6 This artificial numerical diffusion may be further reduced by stronger 
damping of the characteristic radiative wavespeed in the optically thick limit 
(Eq. 1471 . However, it might cause problems in the intermediate regime. 
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Figure 9. The radiative energy density for the optically thick pulse de- 
scribed in Section 14.41 The green line shows the initial profile while the 
other color show the time evolution. Dotted lines show the exact solution of 
Eg.l68l 

adjusted so as to give constant pressure throughout the domain, 

T = T — . (70) 
Po 

The initial radiative energy density is set to the local thermal 
equilibrium value, and the initial velocities and radiative fluxes are 
zero. We apply periodic boundary conditions at the top and bottom 
and outflow boundary conditions at the right border of the domain. 
At the left border we have the external source of light, which we 
specify with E L = 4crT 4 L , F* = 0.99999£ L , T L = 1007V All other 
quantities are set to match the ambient gas. We evolve the system 
with both the Eddington approximation and Ml closure, assuming 
k= X = P- 

Figure [TO] shows the results at t = 10. By this time, the ini- 
tial radiation wave has passed through the domain and the system 
has reached a stationary state. The upper panel shows the solution 
we obtain with the Eddington approximation. Because this closure 
treats all directions equally, radiation readily diffuses into the re- 
gion behind the blob. As a result, there is no shadow behind the op- 
tically thick blob. The lower panel in Figure [10] shows results with 
the the Ml -closure. In contrast to the case of Eddington closure, 
here the x direction is distinguished because the incoming radia- 
tion at the left boundary moves in this direction. The Ml closure is 
designed to keep flux moving parallel to itself in optically thin re- 
gions for F ss E. As a result, a strong shadow develops behind the 
optically thick blob. This test illustrates one of the key differences 
between the Eddington and Ml closure schemes. 

It is appropriate to mention that the excellent performance of 
the Ml closure scheme in this shadow test problem is partly be- 
cause the setup is particularly favorable. First, we have only a single 
source of radiation. Second, the shadow here is aligned well with 
the grid, which helps to minimize diffusion. For other grid orienta- 
tions, the numerical results would exhibit more diffusion, e.g., see 
Section|4~7l 

The Ml closure assumes that the specific intensity is sym- 
metric with respect to the radiative flux, i.e., only one direction 
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Figure 8. Profiles of the radiative energy density (E) for the optically thin radiative pulse test described in Section [4~4| The top panels show its distribution in 
the xy plane at (from left to right) t = 0, 15 and 35. The orange circle in the first plot denotes the initial width and expands at the speed of light to provide the 
expected pulse front location in the other two plots. The bottom panels show the corresponding profiles measured along y = z = line and the 1 /x 2 dependence 
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Figure 10. Results obtained with the shadow test for a single beam of light. 
The upper panel corresponds to the Eddington approximation and the lower 
panel to the Ml closure scheme. Colors denote radiative energy density, 
contours show density (p = 1, 100, 500) and arrows show the direction and 
magnitude of the reduced radiative flux F' / E. 

is distinguished. It means that this approach is supposed to be less 
efficient when multiple sources of light are involved. To test its per- 
formance for such a setup we imple mented a two-bea m test prob- 
lem similar to the one described in ljiang et al.l J2012f1, We set up 
exactly the same initial conditions for gas and radiation as in the 
two previous tests. This time, however, we set up a reflection sym- 
metry at the lower boundary (y = 0) and we impose an inclined 

7 The key princ iple behind the non-relativistic algorithm described in 
Ijiang et all |20l3) is the use of a "Variable Eddington Tensor" (VET). The 
VET is used to close the radiative equations, relating radiative pressure to 
the local radiative energy density. The VET is computed through a sepa- 
rate radiative transfer solver, which calculates (at each time step) the time- 
independent radiation field (using the fixed fluid background of the previ- 
ous timestep as its input). The authors solve the radiative transfer equations 
along a discrete set of rays, and so their scheme accurately captures all 
shadows that can be resolved by these ray angles. The radiation pressure 
obtained from the VET is then used to evolve the MHD fluid equations. 



(F„ = 0.93E , F = -0.37.Eo) beam on the upper boundary and on 
the part (y > .3) of the left boundary. As a result, the domain is ef- 
fectivelly enlighted by two self-crossing beams of light. We plot the 
result of a numerical simulation in Fig.llllwhere the region of neg- 
ative y-coordinates was plotted by reflecting the y-positive data. In 
the region near the left top and bottom corners, where the beams do 
not overlap, the direction of the flux follows the imposed boundary 
condition. In the region of the overlap the radiative energy density 
increases twice (E = 2E ) while the flux becomes equivalent to the 
superposition of the beam-intrinsic fluxes, i.e., it is purely horizon- 
tal and its x-component equals F* = 2Fi = I.86F0 = 0.93E. The 
clump of optically thick gas is, therefore, effectively illuminated 
by a purely horizontal beam. One could expect it creates a paral- 
lel shadow similar to the one obtained in the single beam problem. 
This is, however, not the case. There is an important difference be- 
tween the beam we imposed on the left boundary in the single beam 
problem and the one which develops in the overlapping region. The 
former had F* » E what implied that the specific intensity was al- 
most a S function in the direction of the flux. The latter, however, 
has F v = 0.93F which means that the implied distribution of the 
specific intensity is only an elongated ellipsoid pointing in the di- 
rection of the flux, i.e., photons moving in other directions than the 
direction of the flux are allowed. This has an effect on the shadow 
produced behind the clump. Instead of sharp parallel edges we get 
regions of the partial shadow (penumbra) resulting from these per- 
pendicular photons allowed by the closure when F x < E. The re- 
gion of the total shadow (umbra) is therefore limited by the edges of 
the penumbra and follows the expected shape (compare Fig. 1 1 in 
|jiangetaljj2012h ) to a good accuracy. A significant difference be- 
tween the exact solution and our numerical one arises in the region 
where the penumbrae overlap. One could expect a uniform, trian- 
gular region of no shadow. The Ml closure, however, produces an 
extra narrow horizontal shadow along the x-axis. 



This test shows limits of the Ml closure approach but at the 
same time stresses the fact that, in principle, it does not limit spe- 
cific intensity to one particular direction (assuming only its symme- 
try with respect to the flux). It performs much better than the Ed- 
dington approximation but in the case of multiple sources of light 
it must be used with caution. 
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4.6 Static atmosphere 

An important aspect of radiation in accretion disks is momentum 
transfer between radiation and gas. The Eddington luminosity limit, 
for instance, arises from this interaction. To validate the treatment 
of gas-radiation momentum exchange in our method, we study a 
static atmosphere which is in equilibrium under the combined ac- 
tion of gravity, gas pressure gradient and radiation force. We con- 
sider a polytropic atmosphere on a spherical object. We take the 
optically thin limit and assume that gas-radiation interactions oc- 
cur only through a scattering coefficient, i.e., k = 0, x = *es (equa- 
tionUUl. 

An analytical solution is easily derived for this model prob- 
lem. Because we assume a polytropic equation of state and set 
k = 0, there is no energy equation, and the radial component of 
the momentum equation (r component of equation [24t is all that 
matters. In the non-relativistic limit (r » 2), assuming stationarity 
(d, = 0) and zero velocity (v' = 0), this equation takes the form 



where 



lgP = 1-/ 
p dr r 2 



(71) 



(72) 



Here F in is the radiative flux imposed as a boundary condition at 
the bottom of the atmosphere, r = r- m , and / gives the ratio of the 
radiative to gravitational (or geometrical) forces; / = 1 corresponds 
to the Eddington limit, where the luminosity is L Edd = 47t/k cs and 
the radiative flux is Fj n = Feaa = l/fes'"^- Since radiative energy 
must be conserved, in the stationary state the flux must satisfy F = 
F- m r 2 Jr 2 (non-relativistic limit). 

Equation M\\ may be solved with the help of the polytropic 
equation of state p = Kp r to give, 



0.4pcirr'. At the innermost radius we set p m = 1(T 15 gcirT 3 (op- 
tically thin atmosphere) and T in = 10 6 K. All the velocities were 
initially zero and the radiative energy density E = F in /0.99. Initial 
values of the gas density and temperature in the domain and in the 
ghost cells were assigned based on the analytical solution. We ran 
four models corresponding to four luminosities: 1CT 10 , 0.1, 0.5 and 
1.0 LEdd- Each model was run up to a time t = 2 x 10 9 M, which 
is sufficient to reach relaxed steady state for these optically thin 
atmospheres. 

Figure [T2l shows the results. The top panel shows the density 
profiles corresponding to the four models. Solid lines are the ana- 
lytical solutions and filled circles correspond to the numerical so- 
lutions. The agreement is very good. The higher the luminosity, 
the flatter is the density profile, indicating the effect of the outward 
force due to radiation. For the particular case of the Eddington lu- 
minosity, the density is perfectly constant, reflecting the fact that 
the gravitational force is exactly balanced by radiation and no pres- 
sure gradient in required. We see that the relaxed numerical solu- 
tion is indistinguishable from the analytical solution, confirming 
that KORAL properly handles gas-radiation momentum exchange. 
The plot of residuals at the bottom of the panel indicates that frac- 
tional deviations in the density are below 10~ 4 . Even this small dis- 
crepancy is at least in part because we are comparing a numerical 
solution from a GR code with an analytical solution derived under 
Newtonian physics. 

The middle panel in Figure [12] shows our results for the ra- 
dial radiative flux. Once again, the models behave very well and 
the agreement with the analytical solution is excellent. Finally, the 
bottom panel shows the residual radial velocities (v r /c). These are 
of the order of 10~ 8 (they should be zero), and appear to be mostly 
driven by slight inconsistencies near the boundaries (possibly again 
because of a slight tension between GR and Newtonian physics). 



P = 



TK 



f 



where 



T -x 1 



(r - iy 



(73) 



(74) 



and pi„ is the assumed density at r = r m . The entropy constant K 
is calculated at the bottom of the atmosphere from the assumed gas 
temperature T m . 

We set up a linear grid of 40 points between r = 10 6 and 
1.4 x 10 6 gravitational radii and we solved the problem using MP5 
reconstruction scheme and the standard Ml closure. We scaled 
all quantities to physical units assuming M = 1M Q and k& = 



4.7 Beam of light near BH 

To test the performance of the code for radiation in strong grav- 
itational field, we study propagation of a beam of light in the 
Schwarzschild metric. We consider three models, in each of which 
a beam of light is emitted in the azimuthal direction at a differ- 
ent radius. We decouple gas and radiation by neglecting absorp- 
tions and scatterings {k = x = 0). We run the models on a two- 
dimensional grid with 30 points distributed uniformly in r between 
r- m and r oM (see Table [4] for values) and 60 points distributed uni- 
formly in azimuthal angle cf> between = and n/2. Initially, we 
assign negligibly small values for all primitive quantities, includ- 
ing the radiation energy density and flux. We use outflow boundary 
conditions on all borders except the region covered by the beam at 
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Figure 12. Results obtained with the static atmosphere test. Numerically 
determined profiles and residuals between the numerical and analytical so- 
lutions are plotted for the density (top panel), radial flux (middle panel) and 
radial velocity (bottom panel, residuals only). Colors denote the Eddington 
ratio of the flux boundary condition F m at the bottom of the atmosphere: 
F in = 10 _10 FEdd (red), O.lFridd (orange), 0.5F Edd (light green) and l.0F w 
(dark green) Circles correspond to the numerical solutions and lines show 
the analytical profiles (equation 17 11 . 



Table 4. Model parameters for the light beam tests 



Model 


team 


''in 


''oul 


1 


3.0 ±0.1 


2.5 


3.5 


2 


6.0 ±0.2 


5.3 


7.5 


3 


16.0 ±0.5 


14.0 


20.5 



the equatorial plane (see the range of n> e a m i n Table|4j, where we set 
the radiation temperature to T^am = 10 10 = 10007o and the flux to 
= 0.99999F. Here T is the initial gas and radiation temperature 
of the ambient medium. We use linear reconstruction with 9=1. 

The panels in Figure Q~3] show the results for the three models 
and geodesies of photons at beam boundaries. Consider the right 
panel, which corresponds to Model 3 (Table [4) with the beam cen- 
tered at rbcam = 16. At such a large radius we do not expect signifi- 
cant bending of photon geodesies and this is indeed the case — the 
beam is only slightly bent towards the BH. We also expect the beam 
to be tighly confined, i.e., it should propagate with a nearly constant 



width, as indicated by the two solid green lines, the true geodesies 
of photons at the boundaries of the beam. However, the numerical 
solution shows significant artificial broadening. This is caused by 
numerical diffusion which is significant whenever the radiative flux 
vector is not aligned with the grid geometry, i.e., the beam is tilted 
with respect to the grid axes. 

The middle panel in Figure[T3]shows Model 2, where the beam 
is centered at the marginally stable orbit: rb ea m = 6. At this radius, 
photon geodesies are significantly deviated by gravity, resulting in 
strong curvature in the beam. The numerical beam follows the cor- 
rect trajectory. Moreover, numerical diffusion is lower in this case 
because the curvature brings the beam into closer alignment with 
the grid. There is in addition some real beam divergence because 
the geodesies at different radii within the beam have different cur- 
vatures (see the solid green lines), but this effect is not very signifi- 
cant. 

Finally, the left panel in Figure [T3] shows Model 1, where the 
center of the beam is exactly at the photon orbit: n, C am = 3. An az- 
imuthally oriented ray at this radius is expected to orbit around the 
BH at a constant r. This is seen clearly in the numerical solution. 
Moreover, since the photon geodesic follows the grid, there is prac- 
tically no numerical diffusion. Indeed, there is less diffusion than 
there should be (compare the numerical beam with the two solid 
green lines). The beam should have some divergence as it propa- 
gates around the BH because photons emitted inside r = 3 curve 
inward and will ultimately fall into the BH, while those emitted 
outside r = 3 curve outward and will move towards infinity. The 
simulated beam does not reproduce this physical broadening very 
well. 



4.8 Radiative spherical accretion 

Our last test problem considers radiative spherical accretion 
onto a non-rotat i ng BH . Thi s problem h a s been studied in the 
past by IVitellol d 19841) and iNobili et all dl99ll) and more re- 
cent ly bylRoedig et all d2012h and iFragile et alj J2012h . We fol- 
low [Fragile et al. 1 d2012l) in the setup of our simulations to facili- 
tate comparison with their results. As in their work, we consider 
Thomson scattering and thermal bremsstrahlung, which give the 
following opacity coefficients, 



X 



l.7xl(r a T-" J -m-/pcm- 
k + 0.4pcm _1 , 



(75) 
(76) 



where p is in g cirT 3 and m p is the mass of the proton. Our numeri- 
cal grid spans from r- m = 2.5 to r out = 2x 10 4 and is resolved by 512 
grid points spaced logarithmically following x = log((r- 2.2)/2) 
where the auxiliary variable x is spaced linearly between values 
corresponding to r- m and r olll . We assume a BH mass of 3M Q . For 
the initial state, we choose the mass accretion rate M (see Table [5] 
for values) and set the density profile accordingly, 



M 



4nr 2 u r 



(77) 



where the radial velocity u r is equal to its free fall value u' 
- V2/r. The gas temperature is given by 



P 

Pout 



(78) 



where T ollt is the temperature at the outer radius and T is the adia- 
batic index. The latter is calculated from the radiation to gas pres- 
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Figure 13. Results for Model 1 (left panel), Model 2 (middle), Model 3 (right), involving light beams propagating near a Schwarzschild BH (see Table|4]for 
model details). The BH is at r = (i.e., x = y = 0). The beams are introduced via a boundary condition on the Jt-axis. The beams initially move vertically, 
i.e., in the azimuthal direction. Color indicates the radiation energy density and arrows show the radiative flux as m easured by a ZAMO . The solid green lines 
indicate true geodescis of photons at the beam boundaries. They were calculated using the ray-tracing code GY0T0 I Vincent et al. 201 1). 



Table 5. Models parameters for radiative spherical accretion tests 



Model 




To[K] 


f _ i>rad 


L/ ^Edd 




E1T6 


1.0 


10 6 


1.2 x 10~ 4 


8.73 X 10~ 


-8 


E10T5 


10.0 


10 s 


1.2 x 10~ 7 


3.26 x 10- 


-6 


E10T6 


10.0 


LO 6 


1.2 x 10~ 4 


6.51 x 10- 


-b 


E10T7 


10.0 


10 7 


1.2 x 10" 1 


1.45 x 10- 


-5 


E100T6 


100.0 


10'' 


1.2 x 10~ 4 


2.00 x 10- 


-4 


Model 


names and parameters after iFrasile et 


al.M2012l) 





sure ratio f p = p ra d/ /? gas of the initial state (Table[5]l, 

; =1+ i/2 +2/p ) 

3 



(79) 



,1 + 2/p, 

The radiative energy density is set to E = 3/ p p gas . 

The numerical simulations are run in one (radial) dimension 
with the MP5 reconstruction scheme, and Ml closure for the ra- 
diation. The primitive quantities at the outer boundary are fixed at 
their initial values, as described above. At the inner boundary we 
apply outflow boundary conditions, with the radial velocity fixed at 
the free-fall value and rest mass density and internal energy extrap- 
olated proportional to r~ 3/2 and r~ 3/2r , respectively. Table[5]lists the 
parameter values we used corresponding to five models. The first 
model, E1T6, is characterized by the lowest mass accretion rate 
and is designed to highlight the ability of our scheme to handle 
optically thin me dia. The other four m odels are identical to simula- 
tions described in iFragile et al.ld2012l) . 

Figure [T4lshows the numerical solutions obtained with K0RAL. 
The top-left panel presents profiles of density, which follow the 
initial profile (equationl77t throughout the simulation. The bottom- 
left panel shows the gas temperature. For all but the coldest model, 
E10T5, the temperature follows equation ( 178) . In the case of model 
E10T5, the gas is hotter than the analytical result. This is because 
of gas-radiation coupling which heats up the gas as it approaches 
the BH (the analytical solution assumes that there is no interaction). 
The small kinks in the temperature profiles near the inner bound- 
ary are an artifact of the inner boundary condition. They do not 
influence the rest of the solution since information cannot travel 
upstream in the supersonic flow. 



The top-right and bottom-right panels in Figure[T4]show radial 
profiles of the radiative energy density and flux for the five models. 
Both quantities follow roughly an r~ 2 scaling, reflecting the fact 
that in steady-state (barring redshift factors) the luminosity is equal 
to AnFr 1 and should be conserved. 

Because the flux in these models is non-negligible compared 
to the energy density (e.g., F as 0.9E for the E10 family of models), 
the Eddington closure scheme does not work very well, especially 
at low optical depth. For instance, IFragile et al.l d2012l) used Ed- 
dington closure and obtained unphysical noise or breaks in their 
profiles of radiative quantities (see their Figure 5) in all models 
with M < 300MEdd- This just reflects the fact that their closure 
cannot handle optically thin media. Our algorithm uses the Ml 
closure scheme and has no problems with either optically thick 
or thin regimes. To emphasize this point, we have solved an ad- 
ditional model, E1T6, in which the accretion rate i s an order of 
magni tude lower than the smallest rate considered by Frag ile et al.l 
(2012). K0RAL works fine for this model, and can, in fact, handle 
even more extreme configurations, both at lower and higher accre- 
tion rates. 

For direct c omparison of our results with those reported in 
IFragile eta0 ( l2012l) , we have calculated for all our models the lu- 
minosities, 



L = AnFr 



(80) 



emerging at radius r = 10 3 . Our r esults, shown in Fig ure [T5l are 
consistent with those obtained by IFragile et al.l d2012l) (compare 
their Fi gure 6). Note, howev er, an important difference. The scheme 
used bv lFragile et alj d2012l) is explicit, therefore it cannot reliably 
treat optically-thin media. This causes the luminosities they quote 
to be sensitive to the radius at which they are measured. For in- 
stance, reading off the luminosity values in their simulations (Fig- 
ure 5 of their paper) at a distance at which the flow is optically thin, 
e.g., r = 10 4 , one obtains erroneous lower values of the luminosity. 



5 SUMMARY 

In this paper we have introduced a semi-implicit numerical scheme 
for general relativistic radiation hydrodynamics. The scheme is 
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Figure 14. Numerical results obtained with KORAL for five models of spherical Bondi accretion with radiation. The panels show density (top-left), radiative 
energy density (top-right), gas temperature (bottom-left) and radiative radial flux (bottom-right). Parameters of the models are given in Table|5] The results 
show that the code handles optically-thin and optically-thick regions equally well, without producing unphysical oscillations. 



based on a covariant formulation of the Ml closure scheme for the 
radiation moments. The radiative source terms are handled semi- 
implicitly, and hence this approach can handle practically all op- 
tical depths. The algorithm has been implemented and tested in a 
new GRRHD code KORAL. It can be easily incorporated into any 
general relativistic hydrodynamic or magnetohydrodynamic con- 
servative code. 

Our tests indicate that KORAL works well for a variety of phys- 
ical regimes and geometries: optically thick vs optically thin, gas 
dominated vs radiation dominated, flat space vs curved space. Also, 
as expected, we find that Ml closure has some advantages over 
the standard Eddington closure scheme in the case of optically-thin 
media: namely, it accurately propagates light rays and it is able to 
resolve shadows. 

The semi-implicit radiation scheme implemented in KORAL 
does not overwhelmingly slow down the code. Apart from the fact 
that the inclusion of radiation introduces 4 new conserved quanti- 
ties compared to pure hydrodynamics, the only important compu- 
tational steps that affect performance are (i) calculation of the ra- 
diative characteristic wavespeeds (Section [3. 2) . and (ii) numerical 
solution of the system of four non-linear equations that arise in the 
implicit treatment of radiative source terms (Section 13.31 . We ex- 
pect each of these steps could be speeded up with more effort. How- 
ever, even without further improvements, the code performance is 
already sufficiently good that global GRRHD or GRRMHD simula- 
tions of accretion disks with near- or super-Eddington luminosities 
appear feasible. Meanwhile, it would be of great interest to develop 
closure schemes beyond Ml, e.g., by directly evolving the photon 



distribution function and using it to close the radiation moments, 
for conservative GR codes. 
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Figure 15. Luminosity (in Eddington units) emerging from a spherically 
symmetric radiative Bondi accretion flow onto a non-rotating BH, plotted 
as a function of the dimensionless (reduced) accretion rate in units of the 
Eddington accretion rate, M/Mgdcl- Parameters of the models are given in 
Table [5] Low accretion rates, where the accretion flow becomes optically 
thin, are usually problematic for numerical codes. KORAL has no problem 
handling either low or high accretion rates. 
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APPENDIX A: EXPLICIT-IMPLICIT METHOD FOR THE 
RADIATIVE SOURCE OPERATOR 



In Sect. [33] we described a semi-implicit method for applying the 
radiative source terms C in the lab frame. It requires solving a 
four-dimensional system of non-linear equations. Because of that 
fact, the numerical efficiency is limited. Furthermore, the method 
may sometimes fail to produce a solution, e.g., when the initial 
guess is not close enough to the true solution or when the solu- 
tion at an intermediate iteration becomes unphysical such as having 
F > cE. To handle these situations we have developed the approx- 
imate analytical method described in this Appendix, which is both 
robust and failsafe. We use it as a backup solver for the fiducial 
algorithm. Its limitations are discussed below. 

Let us assume that the advective and geometric source terms 
have already been applied as per steps (i)-(ix) in Section |3~T1 The 
only remaining terms are the radiative forces. For instance, equa- 
tions ([25} require us to time-evolve 



d,(R' v ) = -G Y . 



(Al) 



Let us assume to start with that we treat this term explicitly. Then, 
the update of the conserved quantities is given simply by 

AR' V = -G v At. (A2) 

The right hand side can be rewritten as 

- G v At = -gpyG** At = -g^KG" At, (A3) 

where G a is the radiative four-force in the fluid frame. 

In the same explicit spirit, the fluid-frame source terms in 
equations ( 144b may be written as 



AE 

AF j 



-G' At, 
-G j At. 



(A4) 



These updates of E and F-' generally do not correspond to the up- 
dates of the conserved quantities we search for because the spatial 
and temporal derivatives are mixed when moving from one frame 
to another, i.e., the operator splitting is frame-dependent. However, 
comparing the right hand sides of equations ([A3} and JA4L it is 
clear that the change in the conserved quantities is related to the 



updates of the primitive radiative quantities E and F' calculated 
purely in the fluid frame. From equations (IA21 and dA3t it follows 
that 



AR' V = g, v V a A a P , 



(A5) 



where components of the four-vector of updates of primitive ra- 
diative quantities in the fluid frame, A P (calculated according to 
equation lA4t . are given by the following 4-vector, 



[AE, AF'l 



(A6) 



This quantity may be found using either explicit or implicit meth- 
ods. 

The explicit approach is the simplest and fastest in the opti- 
cally thin regime, but it fails in the optically thick regime. To cal- 
culate A P via the explicit approach, we simply set 



A P = [-G' {n) At,-G J M At], 



(A7) 



where the subscript (n) indicates variables evaluated at the begin- 
ning of the current time step. 

The implicit approach works well for all optical depths, and is 
especially important at high optical depths, but it usually involves 
more computations. In principle, one should solve the following set 
of fluid-frame equations for gas internal energy, momenta, radiative 
energy and fluxes, 



u (n+l) ~ U{n) 



ml L , s — ml , 

(n+1) (n) 



E( n +\) — £(„) 



F' 

(n+1) 



F' 



k(E ( „ +1) - 4<rT*„ +l) )At, (A8) 

(A9) 

-K(E tn+l) - AcrTl +l) ), (A10) 

-xK^t, (All) 

where the subscripts (n) and (n+1) denote values at the beginning 
and the end of the current time step, respectively, and k and x are 
computed from the gas properties at time (n + 1). As a simplifica- 
tion, we assume that k and x do not change significantly during a 
single time step, and we set 

k = K(p,T {n) ), (A 12) 

X = x(p,T { n>). (A13) 

Now, equation iAl It becomes decoupled from the others and may 
be solved directly for the updated fluxes, 



F' 

(n+1) 



F' 



(A14) 



1 + X At 

To find the new value of the radiative energy density one has 
to solve equations ( IA8b and ( IAl 0b together. Taking into account the 
ideal gas equation of state, 



'R. 

p = (T - l)u = —pT, 



(A 15) 



where fi is the mean molecular weight, it is straightforward to com- 
bine these two equations into one quartic equation for E ( „ + iy K0RAL 
solves this equation numerically using the Newton method with E {n) 
as the initial guess. However, analytical solvers for quartic equa- 
tions may also be implemented; in fact, a linearized version of the 
quartic is often quite adequate. Once E and F' have been calcu- 
lated, the four-vector of updates, A P , is given by, 



[£(„+!) - E (n) ,F' 



(n+1) 



F J 1 

(") J - 



(A16) 



This update is applied in step (xi) of the algorithm (see Section [3~TT l 
in case the fiducial lab frame implicit method fails. 

The updates of the conserved quantities calculated using this 
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approach diverge from the fiducial semi-analytical results in the 
limit of large velocites u' » 1 and large optical depths per cell 
t » 1. Therefore, in general, one should use this algorithm only 
as a failsafe backup method. However, if one is confident that the 
problem at hand does not involve such conditions, one might con- 
sider using this method as the default, thereby increasing the code 
efficiency. 
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